classdef Legendre < FilterFather
    properties
        NormalizedPoles = [];
    end
    methods
        
        function this = Legendre(guiData)
            this.calculateFilter(guiData);
            minEdge = min([guiData.Fp guiData.Fa])*2*pi;
            maxEdge = max([guiData.Fp guiData.Fa])*2*pi;
            this.calcData([minEdge maxEdge]);
        end
        
        function calculateFilter(obj,guiData)
            Aa = guiData.Aa;
            Ap = guiData.Ap;
            obj.Approximation = 'Legendre';
            obj.FilterType = guiData.filterType;
            Wp = guiData.Fp*2*pi;
            Wa = guiData.Fa*2*pi;
            [Wp Wa] = obj.adjustBand(Wp,Wa);
            obj.DesnormValue = guiData.DenormVal;
            obj.DesnormType = guiData.DenormType;
            fixOrder = guiData.fixOrder;
            MaxOrder = guiData.order;
%             MaxQ = guiData.maxQ;
            obj.normalization(Wp,Wa);
            
            N = obj.LegendreOrder( Ap, Aa);
            
            if MaxOrder == 0
                obj.Order = N;
            else 
                if(N > MaxOrder || fixOrder == true)
                    N = MaxOrder;
                end
                obj.Order = N;
            end
            
            
            [z p k] = obj.LegendreZPK(obj.Order, Ap);
            [z p k] = ZPKadjust(z,p,0);

            obj.K=k;
            obj.NormalizedPoles = p';
            obj.Poles = p;
            obj.Zeros = z;
%             obj.K = real(prod(obj.Poles)/real(prod(obj.Zeros)));
            obj.frecDesnorm(Aa);
            obj.K = abs(real(prod(obj.Poles)/real(prod(obj.Zeros))));
            obj.desnormalization();
            obj.TransferFunction = obj.TransferFunction*10^(guiData.gain/20);
        end
             
        function n = LegendreOrder(obj, Ap, Aa)

            nMax = 33;
            n = 1;
            while(~obj.CumplePlantilla(n, Ap, Aa) && n < nMax)
            n = n + 1;
            end
    
            if (n == nMax)
                errordlg('Maximum order for Legendre approximation reached. (n = 33)'); 
                n = nMax;
            end
        end
        
        function cumple = CumplePlantilla(obj, n, Ap, Aa)

            [z p k] = obj.LegendreZPK(n, Ap);
            [z p k] = ZPKadjust(z,p,0);
            HL = zpk(z,p,k);
            Magnitud = bode(HL,obj.WaNormalized);
            AaConNactual = -20*log10(Magnitud);
    
    
            if(AaConNactual > Aa)
                cumple = 1;
            else
                cumple = 0;
            end
        end
        
        function [z p k] = LegendreZPK(obj, n, Ap)
            
            if ~rem(n,2) % Si 'n' es par

                k = n/2 - 1;
                b = zeros(1,k+1);
                b0 = 1/sqrt((k+1)*(k+2));
                if ~rem(k,2)
                    for i = 2:2:k
                        b(i) = (2*i+1)*b0;
                    end
                else
                    for i = 1:2:k
                        b(i) = (2*i+1)*b0;
                    end
                    b0 = 0;
                end

                L = b0;
                for i = 1:k
                    Pi = obj.LegendrePoly(i);
                    L = obj.sumarPoly(L, b(i) * Pi);
                end

                L = conv(L,L); %L^2    
                L = conv(L, [1 1]); %(x+1) * L^2

            else

                k = (n-1)/2;
                a0 = 1/(sqrt(2)*(k+1));
                L = a0;
                for i = 1:k
                    Pi = obj.LegendrePoly(i);
                    L = obj.sumarPoly(L, (2*i+1)*a0 * Pi);
                end

                L = conv(L,L);  %L^2

            end

            %Integra
            L = polyint(L);
     
            %%%%%%%%%%%%%% Barrow %%%%%%%%%%%%%%%%%%%
            % limite superior
            or = size(L,2) - 1;
            subsPoly = [2 0 -1];
            Ln = [0];
            for i = or:-1:1
                subsNpoly = subsPoly;
                for j = 1:(i-1)
                    subsNpoly = conv(subsPoly, subsNpoly);    %subs^i
                end
                Ln = obj.sumarPoly(subsNpoly * L(1,or-i+1), Ln);
            end
            Ln = obj.sumarPoly(Ln, L(or+1));

            % limite inferior
            li = zeros(size(Ln));
            li(end) = polyval(L,-1);

            Ln = obj.sumarPoly(Ln, -li);
            Ln((Ln<1e-6) & (Ln>-1e-6)) = 0;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            epsilon = sqrt((10^(Ap/10)-1));

            denominador = obj.sumarPoly(1, epsilon^2 * Ln);

            p = 1j*roots(denominador);
            p = p(real(p)<0);   
            z = [];
            k = abs(prod(p));
        end
        
        function Pn = LegendrePoly(obj,n)

        %Calcula el polinomio de Legendre de N-esimo orden 
        % Rodrigues' formula:
        %                     Pn(x) = 1/(2^n * n!) * d^n/dx^n [(x^2-1)^n]

            %(x^2-1)^n
            polinomioBase = [1 0 -1];
            poln = polinomioBase;
            for i=1:n-1
                poln = conv(polinomioBase, poln); % la convolucion discreta 
            end                                   % calcula el producto de 
                                                  % los polinomios
    
            %d^n/dx^n
            for i=1:n
                poln = polyder(poln);
            end
    
            Pn = 1/(2^n * prod(1:n)) * poln;
        end
       
        function suma = sumarPoly(obj,A,B)

            if length(A) >= length(B)
                N = length(B);
                sumando = B; 
                suma = A;
            else
                N = length(A); 
                sumando = A;
                suma = B; 
            end
   
            suma(end+1-N:end) = suma(end+1-N:end)+sumando;
        end
        
    end

end
